Image processing of video using noise estimate

ABSTRACT

A method of processing image date representing an image of a scene to generate an estimate of noise present in the image data. The method comprises evaluating a function for different values of the estimate, the function taking as input an estimate of the noise, and determining an estimate of the noise for which the function has an optimum value.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application is a divisional of U.S. Ser. No. 11/908,736, filed Sep.14, 2007, now U.S. Pat. No. ______, which is a §371 National Stagefiling of PCT/GB2006/000974, filed Mar. 16, 2006, and which claims thebenefit of U.S. Ser. No. 60/665,491, filed Mar. 25, 2005. The entirecontents of these prior applications are hereby incorporated byreference.

TECHNICAL FIELD

The present invention relates to image processing methods, and moreparticularly but not exclusively to methods for enhancing images whichare at least partially affected by the effects of noise.

BACKGROUND OF THE INVENTION

There are many sources of noise which may degrade an image of a scene.For example an image of a scene will often be degraded by opticalscattering of light caused, for example by fog or mist. This opticalscattering results in additional lightness being present in some partsof the image, and has been referred to as “airlight” in the literature.It is desirable to process an image so as to remove components of pixelvalues which are attributable to airlight.

If the distance between a camera position and all points of a scenerepresented by an image generated by the camera is approximatelyconstant, airlight can be estimated and removed by applying equation (1)to each pixel of the image:

y=m(x−c)  (1)

where:

-   -   x is an original pixel value;    -   c is a correction selected to represent “airlight”;    -   m is a scaling parameter; and    -   y is a modified pixel value.

Assuming that parameter c is correctly chosen, processing each pixel ofa monochrome image in accordance with equation (1) will enhance an imageby removing air light. However determination of an appropriate value forthe parameter c is often problematic.

Various known methods exist for estimating the parameter c by usingcontrast measurements such as a ratio of the standard deviation of pixelvalues to the mean of pixel values. However, such contrast measures donot discriminate between airlight-induced contrast loss and inherentlylow contrast scenes. For example an image of sand dunes would oftenprovide little contrast between the light and dark parts of the scene,even when no airlight is present. Thus ad-hoc schemes to determine theparameter c will sometimes result in severe image distortion.

The method described above with reference to equation (1) is applicableto a monochrome image. Further problems arise if a colour image is to beprocessed. Typically the airlight contribution to a pixel value, andhence the value of the parameter c, will depend upon the wavelength(colour) of the light. Thus, if equation (1) is to be applied to colourimages, different values of the parameter c may be needed for red, greenand blue channels of the image.

The methods described above assume that the camera position isequidistant from all points in a scene represented by an image.Published European Patent EP 0839361 describes a method developed by oneof the present inventors, in which different values of the parameter cin equation (1) are used for different pixels, in dependence upondistance between the camera position and a position in a scenerepresented by that pixel. This invention arose from a realisation thatbackscattered light may vary in dependence upon the distance between acamera position and a position in a scene.

SUMMARY OF THE INVENTION

It is an object of embodiments of the present invention to obviate ormitigate at least some of the problems outlined above.

According to the present invention, there is provided, a method ofprocessing image data representing an image of a scene to generate anestimate of noise present in the image data. The method involves afunction taking as input an estimate of the noise. The function isevaluated for a plurality of different values of the estimate, and anestimate for which the function has an optimum value is determined. Theoptimum value may be a stationary point of a function, for example aminimum.

Thus, the present invention provides a robust method of determiningnoise present in an image by simply evaluating a function anddetermining a value of the noise estimate for which the function has anoptimum value.

The term noise as used in this document is used to mean any artefact ofan image which imparts slow varying offsets to that image. Examples ofnoise are described in further detail below although it will beappreciated that examples of such noise include airlight.

The image data may comprise pixel values for each of a plurality ofpixels of the image, and the function may take as input at least asubset of the pixel values.

The method may further comprise filtering the image data to generatefiltered image data comprising filtered pixel values for each of theimage pixels. The function may then take as input at least a subset ofsaid filtered pixel values. The subset is preferably selected so as tobe the same subset of pixels for which the function takes as input pixelvalues.

The function described above may be of the form

$\begin{matrix}{{S(\lambda)} = {\frac{1}{K}{\sum\limits_{k = 1}^{k = K}{{\left( \frac{p_{k} - \overset{\_}{p_{k}}}{\overset{\_}{p_{k}} - \lambda} \right)^{2} \cdot \exp}\frac{1}{K}{\sum\limits_{k = 1}^{k = K}\left( {\ln \left( {\overset{\_}{p_{k}} - \lambda} \right)}^{2} \right)}}}}} & (2)\end{matrix}$

where:

-   -   K is the number of pixels to be processed;    -   p_(k) is the value of pixel k    -   p_(k) is the value of pixel k after application of the lowpass        filter described above;    -   λ is the value of the parameter; and    -   S(λ) is the function to be optimised.

The optimum value may be determined using any of the large number ofnumerical analysis techniques which are well known in the art. Suchnumerical analysis techniques preferably begin with an initial estimateof zero for the noise included in the image.

It is known that some noise present within images varies in dependenceon the distance between a camera position and a position in the scene.Such variants may be taken into account using the methods set out aboveby generating a plurality of different estimates.

The invention also provides a method of removing noise from an image.The method comprises estimating noise included in the image using amethod as set out above, and then processing each pixel of the image toremove the estimate of said noise to generate said output pixel values.

The method of removing noise may further comprise multiplying the outputpixel values by a predetermined coefficient. Multiplying output pixelvalues by a predetermined coefficient may provide illuminationcompensation using one of a number of well known illuminationcompensation techniques.

The noise included in the image maybe at least partially attributable toone or more of atmospheric back scattered light, a dark current effectof a camera, or dirt on a camera lens.

According to a further aspect of the present invention, there isprovided, a carrier medium carrying computer readable instructionsconfigured to cause a computer to carry out a method according to anyproceeding claims.

The invention further provides a computer apparatus for generating anestimate of noise present in image data representing an image. Theapparatus comprises a program memory containing processor readableinstructions, and a processor configured to read and executeinstructions stored in said program memory. The processor readableinstructions comprise instructions configured to cause the computer tocarry out a method as described above.

According to a further aspect of the present invention there isprovided, a method of processing a plurality of frames of video data togenerate enhanced video data. The method comprises processing a firstframe of video data to generate an estimate of noise included withinsaid frame of video data, storing data indicative of said estimate ofnoise, processing at least one further frame of video data to removenoise from said further frame of video data using said stored estimateof noise included in said first frame of video data and outputting anenhanced frame of video data generated using said further frame of videodata.

Thus, by using the methods set out above it is possible for processingof a first frame of video data to proceed relatively slowly whileprocessing of the at least one further frame of video data may takeplace relatively quickly so as to prevent artefacts of ‘jerkiness’ beingapparent in the output video data.

The first frame of video data may be processed in accordance with themethods set out above. This processing may further generate a contrastenhancement parameter.

BRIEF DESCRIPTION OF THE DRAWINGS

Embodiments of the present invention will now be described, by way ofexample, with reference to the accompanying drawings, in which:

FIG. 1 is a flowchart of an image processing method in accordance withthe present invention;

FIG. 2 is an example image to which processing in accordance with FIG. 1is applied;

FIG. 3 is a graph showing values of a parametric function used in theprocessing of FIG. 1 when applied to the image of FIG. 2;

FIG. 4 is an illustration of scene geometry; and

FIG. 5 is a schematic illustration of an implementation of an imageprocessing method in accordance with the invention.

DETAILED DESCRIPTION OF THE ILLUSTRATED EMBODIMENTS

Pixel values of an image may include a component attributable to noisecaused by atmospheric backscattered light, often referred to asairlight. Such noise can be removed from the image using equation (1)set out above, and recalled below:

y=m(x−c)  (1)

where:

-   -   x is an original pixel value;    -   c is a correction selected to represent “airtight”;    -   m is a scaling parameter; and    -   y is a modified pixel value.

An embodiment of the present invention provides a method for estimatinga value of the parameter c. This method is illustrated in the flowchartof FIG. 1. The method operates by sampling pixel values associated withat least a subset of the pixels defining the image, inputting thesevalues into a parametric function, and determining a value for theparameter of the parametric function for which the function has anoptimum value. The value of the parameter so determined is the valuegenerated for the parameter c in equation (1).

Referring to FIG. 1, at step S1 pixel values p associated with at leastsome of the pixels within the image to be processed are sampled. At stepS2 a low-pass filter is applied to the image to be processed to generatefiltered pixel values p for each of the image pixels. The filter ispreferably a low-pass filter having a Gaussian shaped kernel and a valueof sigma equal to approximately five. The generation of such a low passfilter will be well known to one of ordinary skill in the art.

At step S3 the parameter λ of the parametric function (described below)is initialised, typically to a value of 0. The parametric functionevaluated at step S4 of FIG. 1 is that of equation (2):

$\begin{matrix}{{S(\lambda)} = {\frac{1}{K}{\sum\limits_{k = 1}^{k = K}{{\left( \frac{p_{k} - \overset{\_}{p_{k}}}{\overset{\_}{p_{k}} - \lambda} \right)^{2} \cdot \exp}\frac{1}{K}{\sum\limits_{k = 1}^{k = K}\left( {\ln \left( {\overset{\_}{p_{k}} - \lambda} \right)}^{2} \right)}}}}} & (2)\end{matrix}$

where:

-   -   K is the number of pixels to be processed;    -   p_(k) is the value of pixel k    -   p_(k) is the value of pixel k after application of the lowpass        filter described above;    -   λ is the value of the parameter; and    -   S(λ) is the function to be optimised.

As mentioned above, λ is initialised to a value of 0 at step S3 andequation (2) is then evaluated at step S4. At step S5 a check is made todetermine whether the evaluated value of equation (2) is an optimumvalue, in this case a minimum value. When it is determined that this isthe case, parameter c is set to be the determined value of λ at step S6.Until a value of λ for which S(λ) is a minimum is found, λ is updated atstep S7, and processing returns to step S4.

Location of the minimum value of S(λ) (step S3 to S7) can be carried outusing a number of known numerical optimisation techniques. Furtherdetails of one such technique are presented below.

It has been indicated that processing is carried out to determine avalue of λ, for which the function S(λ) is a minimum. It is now shownthat such a value of λ will be a good estimate for the parameter c inequation (1).

Any image can be regarded as a collection of regions with each regioncorresponding to some kind of visible surface. For example a tree may beseen as collection of foliage, trunk and branch regions. Each region hasa certain spectral reflectivity which is associated with its localpigmentation. Within each region there is a variation in brightness thatis mainly due to macroscopic variations in local surface orientation.This consideration of an image corresponds to a fairly general and basicmodel for image texture.

A key feature of the model described above is that the fractionalvariation in brightness is approximately independent of the actualbrightness. Any offsets due to airlight will change this feature and soare detectable. This is used as a basis for the processing describedabove with reference to FIG. 1 and equation (2).

In the described embodiment the local texture variation is assumed tohave a Normal (Gaussian) distribution.

An image to be processed is assumed to contain M regions and the meanbrightness of the a region m is R_(m), where m=1 . . . M. No loss ofgenerality is implied as M could be very large and the regions verysmall. An explicit segmentation of the image into these M regions is notperformed.

The total number of pixels in the image is denoted by K, and thefraction of those pixels in region m is denoted by K_(m). The k_(th)pixel in region m is denoted by p_(k) ^(m).

The k_(th) pixel in region m (p_(k) ^(m)) is represented as follows:

p _(k) ^(m) =R _(m)(1+N)+c,  (3)

where:

-   -   N=N(0,σ) is a Normal random variable with zero mean, standard        deviation σ;    -   c is a scalar constant indicative of airlight.

The low-pass filter described above smoothes out random fluctuations inthe local pixel value p_(k) ^(m) to produce a “smoothed” value p_(k)^(m) for each pixel.

If the spatial extent of the filter is sufficiently large for thepurpose of smoothing but still small with respect to the size of theimage region (so that border effects are not significant) then:

p _(k) ^(m) ≈R _(m) +c.  (4)

That is, the filtering removes the variation represented by normalvariable N.

Recalling equation (2):

$\begin{matrix}{{S(\lambda)} = {\frac{1}{k}{\sum\limits_{k = 0}^{k < K}{{\left( \frac{p_{k} - \overset{\_}{p_{k}}}{\overset{\_}{p_{k}} - \lambda} \right)^{2} \cdot \exp}\frac{1}{k}{\sum\limits_{k = 0}^{k < K}\left( {\ln \left( {\overset{\_}{p_{k}} - \lambda} \right)}^{2} \right)}}}}} & (2)\end{matrix}$

While equation (2) requires summations over all pixels of the image, itcan be rewritten in terms of each region m of the M regions, recallingthat K_(m), is the fraction of pixels in region m.

Recalling also equations (3) and (4), equation (2) can be written as:

$\begin{matrix}{{S(\lambda)} = {\sum\limits_{m = 1}^{m = M}{{K_{m}\left( \frac{\left( {\left( {{R_{m}\left( {1 + N} \right)} + c} \right) - \left( {R_{m} + c} \right)} \right)^{2}}{\left( {R_{m} + c - \lambda} \right)^{2}} \right)} \cdot {\exp \left( {\sum\limits_{m = 1}^{m = M}\left( {K_{m} \cdot {\ln \left( {R_{m} + c - \lambda} \right)}^{2}} \right)} \right)}}}} & (5)\end{matrix}$

Simplifying equation (5), and recalling that N=N(0,σ) gives:

$\begin{matrix}{{S(\lambda)} = {\sum\limits_{m = 1}^{m = M}{{K_{m}\left( \frac{R_{m}^{2}\sigma^{2}}{\left( {R_{m} + c - \lambda} \right)^{2}} \right)} \cdot {\exp \left( {\sum\limits_{m = 1}^{m = M}\left( {K_{m} \cdot {\ln \left( {R_{m} + c - \lambda} \right)}^{2}} \right)} \right)}}}} & (6)\end{matrix}$

given that σ² is the expectation value of the square of N.

Define:

u=c−λ.  (7)

and substitute equation (7) into equation (6):

$\begin{matrix}{{S\left( {c - u} \right)} = {{\sigma^{2}\left( {\sum\limits_{m = 1}^{m = M}{K_{m}\frac{R_{m}^{2}}{\left( {R_{m} + u} \right)^{2}}}} \right)} \cdot {\exp \left( {2{\sum\limits_{m = 1}^{m = M}{K_{m}{\ln \left( {R_{m} + u} \right)}}}} \right)}}} & (8)\end{matrix}$

Define:

$\begin{matrix}{f = {\sum\limits_{m = 1}^{m = M}{K_{m}\frac{R_{m}^{2}}{\left( {R_{m} + u} \right)^{2}}}}} & (9) \\{g = {\exp \left( {2{\sum\limits_{m = 1}^{m = M}{K_{m}{\ln \left( {R_{m} + u} \right)}}}} \right)}} & (10)\end{matrix}$

Differentiating equation (9) with respect to u gives:

$\begin{matrix}{\frac{f}{u} = {{- 2}{\sum\limits_{m = 1}^{m = M}{K_{m}\frac{R_{m}^{2}}{\left( {R_{m} + u} \right)^{3}}}}}} & (11)\end{matrix}$

Differentiating equation (10) with respect to u gives:

$\begin{matrix}{\frac{g}{u} = {2{{\exp \left( {2{\sum\limits_{m = 1}^{m = M}{K_{m}{\ln \left( {R_{m} + u} \right)}}}} \right)} \cdot {\sum\limits_{m = 1}^{m = M}{K_{m}\frac{1}{R_{m} + u}}}}}} & (12)\end{matrix}$

From equation (8):

S(λ)=σ² f·g  (13)

Therefore:

$\begin{matrix}{\frac{S}{u} = {{\sigma^{2}{\left( {f \cdot g} \right)}} = {\sigma^{2}\left( {{{f} \cdot g} + {f \cdot {g}}} \right)}}} & (14)\end{matrix}$

Substituting equations (9), (10), (11) and (12) into equation (14)gives:

$\begin{matrix}{\frac{S}{u} = {{{- 2}{\sigma^{2}\left( {\sum\limits_{m = 1}^{m = M}{K_{m}\frac{R_{m}^{2}}{\left( {R_{m} + u} \right)^{3}}}} \right)}{\exp \left( {2{\sum\limits_{m = 1}^{m = M}{K_{m}{\ln \left( {R_{m} + u} \right)}}}} \right)}} + {2{\sigma^{2}\left( {\sum\limits_{m = 1}^{m = M}{K_{m}\frac{1}{R_{m} + u}}} \right)}{\exp \left( {2{\sum\limits_{m = 1}^{m = 1}{K_{m}{\ln \left( {R_{m} + u} \right)}}}} \right)}\left( {\sum\limits_{m = 1}^{m = M}{K_{m}\frac{R_{m}^{2}}{\left( {R_{m} + u} \right)^{2}}}} \right)}}} & (15)\end{matrix}$

Differentiating equation (7) with respect to λ gives:

$\begin{matrix}{\frac{u}{\lambda} = {- 1}} & (16)\end{matrix}$

It can be stated that:

$\begin{matrix}{\frac{S}{\lambda} = {\frac{S}{u} \cdot \frac{u}{\lambda}}} & (17)\end{matrix}$

That is

$\begin{matrix}{\frac{S}{\lambda} = {{2{\sigma^{2}\left( {\sum\limits_{m = 1}^{m = M}{K_{m}\frac{R_{m}^{2}}{\left( {R_{m} + u} \right)^{3}}}} \right)}{\exp \left( {2{\sum\limits_{m = 1}^{m = M}{K_{m}{\ln \left( {R_{m} + u} \right)}}}} \right)}} - {2{\sigma^{2}\left( {\sum\limits_{m = 1}^{m = M}{K_{m}\frac{1}{R_{m} + u}}} \right)}{\exp \left( {2{\sum\limits_{m = 1}^{m = 1}{K_{m}{\ln \left( {R_{m} + u} \right)}}}} \right)}\left( {\sum\limits_{m = 1}^{m = M}{K_{m}\frac{R_{m}^{2}}{\left( {R_{m} + u} \right)^{2}}}} \right)}}} & (18)\end{matrix}$

It has been stated above that λ=c when S is a minimum, mat is when

$\frac{S}{\lambda} = 0.$

If λ=c, then from equation (7) u=0. Substituting u=0 into equation (18)gives:

$\begin{matrix}{\frac{S}{\lambda} = {{2{\sigma^{2}\left( {\sum\limits_{m = 1}^{m = M}{K_{m}\frac{1}{R_{m}}}} \right)}{\exp \left( {2{\sum\limits_{m = 1}^{m = M}{K_{m}\ln \; R_{m}}}} \right)}} - {2{\sigma^{2}\left( {\sum\limits_{m = 1}^{m = M}{K_{m}\frac{1}{R_{m}}}} \right)}{\exp \left( {2{\sum\limits_{m = 1}^{m = M}{K_{m}\ln \; R_{m}}}} \right)}\left( {\sum\limits_{m = 1}^{m = M}K_{m}} \right)}}} & (19)\end{matrix}$

Given that

${\left( {\sum\limits_{m = 1}^{m = M}K_{m}} \right) = 1},$

it can be seen that

$\frac{S}{\lambda} = 0.$

This supports the claim that equation (2) has a stationary point whenλ=c. This means that c may be determined by applying an efficientnumerical search procedure to the equation (2).

Having shown that a value of λ for which S(λ) is a minimum will providea good estimate for c, a method for minimization of the function S(λ) isnow described. The simplest way to minimize the function of equation (2)is to directly evaluate values of equation (2). However, thisformulation of the cost function suffers from numerical inaccuracies. Ingeneral terms, for numerical optimization algorithms to work efficientlythe target function must be a continuous and smooth function of itsparameters. An alternative, but equivalent, formulation of the costfunction is described below. This formulation is preferred because itleads to more rapid determination of the optimum value of λ. An exampleof a suitable optimization code is the routine E04UCC from the NumericalAlgorithms Group Ltd, Oxford.

Image data is processed on a row by row basis for each row i. Theresults of this row by row processing are then combined in a recursivemanner.

From equation (2):

$\begin{matrix}{{s(\lambda)} = {{\overset{\_}{\left( \frac{P - \overset{\_}{p}}{p - \lambda} \right)^{2}} \cdot \exp}\overset{\_}{\left( {\ln \left( {\overset{\_}{p} - \lambda} \right)}^{2} \right.}}} & (20)\end{matrix}$

Then define:

$\begin{matrix}{{{RS}(i)} = {\frac{1}{M}{\sum\limits_{m = 1}^{M}\frac{\left( {p_{m} - \overset{\_}{p_{m}}} \right)^{2}}{\left( {\overset{\_}{p_{m}} - \lambda_{t}} \right)^{2}}}}} & (21) \\{{{RL}(i)} = {\frac{1}{M}{\sum\limits_{m = 1}^{M}{\ln \left( {\overset{\_}{p_{m}} - \lambda_{t}} \right)}^{2}}}} & (22)\end{matrix}$

for a given row i, where M is the number of columns and m is the columncurrently being processed.

$\begin{matrix}{{{RTL}(n)} = {\frac{1}{n}{\sum\limits_{i = 1}^{n}{{RL}(i)}}}} & (23)\end{matrix}$

Where n is the number of rows. Then it can be seen that:

$\begin{matrix}{{S(n)} = {{\frac{1}{n}\left\lbrack {\sum\limits_{i = 1}^{n}{{RS}(i)}} \right\rbrack} \cdot {\exp \left( {{RTL}(n)} \right)}}} & (24) \\{{S\left( {n + 1} \right)} = {{\frac{1}{n + 1}\left\lbrack {\sum\limits_{i = 1}^{n + 1}{{RS}(i)}} \right\rbrack} \cdot {\exp\left( {{RTL}\left( {n + 1} \right)} \right.}}} & (25)\end{matrix}$

Recalling the definition of RTL(n) from equation (23) equation 25 can berewritten as:

$\begin{matrix}{{S\left( {n + 1} \right)} = {{\frac{1}{n + 1}\left\lbrack {\sum\limits_{i = 1}^{n + 1}{{RS}(i)}} \right\rbrack} \cdot {\exp \left\lbrack {\frac{1}{n + 1}{\sum\limits_{i = 1}^{n + 1}{{RL}(i)}}} \right\rbrack}}} & (26)\end{matrix}$

Rearranging the summations of equation 26 gives:

$\begin{matrix}{{S\left( {n + 1} \right)} = {{{\frac{1}{n + 1}\left\lbrack {\left( {\sum\limits_{i = 1}^{n}{{RS}(i)}} \right) + {{RS}\left( {n + 1} \right)}} \right\rbrack} \cdot \exp}\left\{ {\frac{1}{n + 1}\left\lbrack {\left( {\sum\limits_{i = 1}^{n}{{RL}(i)}} \right) + {{RL}\left( {n + 1} \right)}} \right\rbrack} \right\}}} & (27)\end{matrix}$

Rearranging gives:

$\begin{matrix}{{S\left( {n + 1} \right)} = {\left\lbrack {\left( {\frac{1}{n + 1}{\sum\limits_{i = 1}^{n}{{RS}(i)}}} \right) + \frac{{RS}\left( {n + 1} \right)}{n + 1}} \right\rbrack \cdot {\exp \left\lbrack {\left( {\frac{1}{n + 1}{\sum\limits_{i = 1}^{n}{{RL}(i)}}} \right) + \frac{{RL}\left( {n + 1} \right)}{n + 1}} \right\rbrack}}} & (28)\end{matrix}$

Given that the exponential of a sum is equal to the multiplication ofthe exponentials of the sum's components, equation 28 can be rewrittenas:

$\begin{matrix}{{S\left( {n + 1} \right)} = {\left\lbrack {\left( {\frac{1}{n + 1}{\sum\limits_{i = 1}^{n}{{RS}(i)}}} \right) + \frac{{RS}\left( {n + 1} \right)}{n + 1}} \right\rbrack \cdot {\exp \left\lbrack {\frac{1}{n + 1}{\sum\limits_{i = 1}^{n}{{RL}(i)}}} \right\rbrack} \cdot {\exp \left\lbrack \frac{{RL}\left( {n + 1} \right)}{n + 1} \right\rbrack}}} & (29)\end{matrix}$

Which can be rewritten as:

$\begin{matrix}{{S\left( {n + 1} \right)} = {\left\lbrack {\left( {{\frac{n}{n + 1} \cdot \frac{1}{n}}{\sum\limits_{i = 1}^{n}{{RS}(i)}}} \right) + \frac{{RS}\left( {n + 1} \right)}{n + 1}} \right\rbrack \cdot {\exp \left\lbrack {{\frac{n}{n + 1} \cdot \frac{1}{n}}{\sum\limits_{i = 1}^{n}{{RL}(i)}}} \right\rbrack} \cdot {\exp \left\lbrack \frac{{RL}\left( {n + 1} \right)}{n + 1} \right\rbrack}}} & (30)\end{matrix}$

Which, recalling the definition of RTL(n) from equation (23) gives:

$\begin{matrix}{{S\left( {n + 1} \right)} = {\left\lbrack {\left( {{\frac{n}{n + 1} \cdot \frac{1}{n}}{\sum\limits_{i = 1}^{n}{{RS}(i)}}} \right) + \frac{{RS}\left( {n + 1} \right)}{n + 1}} \right\rbrack \cdot {\exp \left\lbrack {\frac{n}{n + 1}{{RTL}(n)}} \right\rbrack} \cdot {\exp \left\lbrack \frac{{RL}\left( {n + 1} \right)}{n + 1} \right\rbrack}}} & (31)\end{matrix}$

Multiplying out of the expression of equation (31) gives:

$\begin{matrix}{{S\left( {n + 1} \right)} = {{\frac{n}{n +} \cdot {\frac{1}{n}\left\lbrack {\sum\limits_{i = 1}^{n}{{RS}(i)}} \right\rbrack} \cdot {\exp \left\lbrack {\frac{n}{n + 1}{{RTL}(n)}} \right\rbrack} \cdot {\exp \left\lbrack \frac{{RL}\left( {n + 1} \right)}{n + 1} \right\rbrack}} + {\frac{{RS}\left( {n + 1} \right)}{n + 1} \cdot {\exp \left\lbrack \frac{{RL}\left( {n + 1} \right)}{n + 1} \right\rbrack} \cdot {\exp \left\lbrack {\frac{n}{n + 1}{{RTL}(n)}} \right\rbrack}}}} & (32)\end{matrix}$

Which can be rewritten as:

$\begin{matrix}{{S\left( {n + 1} \right)} = {{\frac{n}{n + 1} \cdot {\frac{1}{n}\left\lbrack {\sum\limits_{i = 1}^{n}{{RS}(i)}} \right\rbrack} \cdot {\exp \left\lbrack {\left( {1 + \frac{- 1}{n + 1}} \right) \cdot {{RTL}(n)}} \right\rbrack} \cdot {\exp \left\lbrack \frac{{RL}\left( {n + 1} \right)}{n + 1} \right\rbrack}} + {\frac{{RS}\left( {n + 1} \right)}{n + 1} \cdot {\exp \left\lbrack \frac{{RL}\left( {n + 1} \right)}{n + 1} \right\rbrack} \cdot {\exp \left\lbrack {\frac{n}{n + 1}{{RTL}(n)}} \right\rbrack}}}} & (33)\end{matrix}$

Which can be rewritten as:

$\begin{matrix}{{S\left( {n + 1} \right)} = {{\frac{n}{n + 1} \cdot {\frac{1}{n}\left\lbrack {\sum\limits_{i = 1}^{n}{{RS}(i)}} \right\rbrack} \cdot {\exp \left\lbrack {{{RTL}(n)} + {\frac{- 1}{n + 1}{{RTL}(n)}}} \right\rbrack} \cdot {\exp \left\lbrack \frac{{RL}\left( {n + 1} \right)}{n + 1} \right\rbrack}} + {\frac{{RS}\left( {n + 1} \right)}{n + 1} \cdot {\exp \left\lbrack \frac{{RL}\left( {n + 1} \right)}{n + 1} \right\rbrack} \cdot {\exp \left\lbrack {\frac{n}{n + 1}{{RTL}(n)}} \right\rbrack}}}} & (34)\end{matrix}$

Rewriting the second term of the multiplication gives:

$\begin{matrix}{{S\left( {n + 1} \right)} = {{\frac{n}{n + 1} \cdot {\frac{1}{n}\left\lbrack {\sum\limits_{i = 1}^{n}{{RS}(i)}} \right\rbrack} \cdot {\exp \left\lbrack {{RTL}(n)} \right\rbrack} \cdot {\exp \left\lbrack {\frac{- 1}{\; {n + 1}}{{RTL}(n)}} \right\rbrack} \cdot {\exp \left\lbrack \frac{{RL}\left( {n + 1} \right)}{n + 1} \right\rbrack}} + {\frac{{RS}\left( {n + 1} \right)}{n + 1} \cdot {\exp \left\lbrack \frac{{RL}\left( {n + 1} \right)}{n + 1} \right\rbrack} \cdot {\exp \left\lbrack {\frac{n}{n + 1}{{RTL}(n)}} \right\rbrack}}}} & (35)\end{matrix}$

Recalling the definition of S(n) from equation (24) gives:

$\begin{matrix}{{S\left( {n + 1} \right)} = {{\frac{n}{n + 1} \cdot {S(n)} \cdot {\exp \left\lbrack {\frac{- 1}{n + 1}{{RTL}(n)}} \right\rbrack} \cdot {\exp \left\lbrack \frac{{RL}\left( {n + 1} \right)}{n + 1} \right\rbrack}} + {\frac{{RS}\left( {n + 1} \right)}{n + 1} \cdot {\exp \left\lbrack \frac{{RL}\left( {n + 1} \right)}{n + 1} \right\rbrack} \cdot {\exp \left\lbrack {\frac{n}{n + 1}{{RTL}(n)}} \right\rbrack}}}} & (36)\end{matrix}$

Thus, it can be seen that the cost function of equation 2 can beeffectively evaluated using equation (36). That is, by suitableinitialization of the variables RS(i) and RL(i) and RTL(n) for the firstrow of the image the cost function S can be calculated on a row by rowbasis as indicated by equation 36.

In a computer program implementation of the algorithm a for loop processeach row of the image in turn. An initial value of RS(i) is calculatedfor the first row and stored in a variable RowProdSum. An initial valuefor RL(i) for the first row is also calculated and stored in a variableRow_Sum_Log. A variable S is then initialized and updated for each rowwith equation 23 being used to update values of RTL as appropriate. Avariable wf is set to be equal 1.0/(i+1.0) which corresponds to the term

$\frac{1}{n + 1}$

in the above equations. Appropriate program code is shown in codefragment 1 below.

Two double precision floating point variables, A and RTB (whichrepresents RTL), are used in the line-recursive structure shown below.The line number is denoted by i.

CODE FRAGMENT 1   {Process each image row, starting at row 0. The rownumber is i} if (i==0) {  RTB = row_sum_log;  double exp_term =exp(RTB);  A = RowProdSum*exp_term; } else {  double wf = 1.0/(i+1.0); double exp_term = exp(wf*row_sum_log);  double LineProdSumProd =wf*RowProdSum*exp_term;  double RTA1_term1 = wf*i*A*exp(−wf*RTB); double RTA1_term2 = exp(wf*i*RTB);  A = RTAl_term1*exp_term +LineProdSumProd*RTA1_term2;  RTB = wf*(RTB*i + row_sum_log [colour]); }i = i + 1; {when all rows are processed the final cost function value isA}

Application of the method illustrated in FIG. 1 to an example imageshown as FIG. 2 is now described.

The image of FIG. 2 is a synthetic image that has been formed toillustrate the method of the described embodiment of the invention. Theimage comprises two regions with different intensity levels.

Initially an image was generated having an intensity of 30 in an area 1and intensity of 120 in an area 2. A texture effect was simulated byadding a random number to each pixel. The random number was generatedusing a Normal distribution with zero mean and a standard deviationequal to one-tenth of the pixel value. The effect of airtight wassimulated by adding a constant value of 50 to each pixel. In the imageof FIG. 2 the mean value in the area 1 is 80 and the mean value in thearea is 150.

FIG. 3 shows a plot of the function S(λ) defined in equation (2) fordiffering values of λ The values of S(λ) are generated using image datataken from all pixels of FIG. 2, and pixel values after application ofan appropriate lowpass filter in the manner described above. It can beseen that S(λ) has a minimum value when λ=53.42. This is close to theaccurate value of λ=50, (50 being the simulated airtight). Slightdistortion along the edge of two patches results in the estimated offsetbeing slightly offset from the true value of λ=50.

When pixel values in the vicinity of the edge are excluded fromevaluation of S(λ), the function has a minimum at λ=49.9, which is veryclose to the true value. Elimination from the calculation of pixels neara sharp edge in an image can be carried out automatically by computingan edge strength measure (for example the modulus of the gradient) andexcluding strong edges from the analysis.

The method set out above will work for any image in which there are atleast two regions of different intensity. This is not restrictive inpractice.

Equation (2) set out above can be modified for situations where theimage contains a known level of random noise. In such a case S(λ) isdefined by equation (37) below:

$\begin{matrix}{{S(\lambda)} = {\frac{1}{K}{\sum\limits_{k = 1}^{k = K}{\left( \frac{\left( {p_{k} - \overset{\_}{p_{k}}} \right)^{2}}{\left( {\overset{\_}{p_{k}} - \lambda} \right)^{2} + \frac{A_{v}}{T_{v}}} \right)\exp \frac{1}{K}{\sum\limits_{k = 1}^{k = K}\left( {{\ln \left( {\overset{\_}{p_{k}} - \lambda} \right)}^{2} + \frac{A_{v}}{T_{v}}} \right)}}}}} & (37)\end{matrix}$

where:

-   -   K is the number of pixels to be processed;    -   p_(k) is the value of pixel k    -   p_(k) is the value of pixel k after application of the lowpass        filter described above;    -   A_(v) denotes the variance of the additive noise component;    -   T_(v) is a constant related to the image texture;    -   λ is the value of the parameter to be determined; and    -   S(λ) is the function to be optimised.

Again, the value of λ for which S(λ) is a minimum will be a goodestimate of the airlight component present in the image.

The preceding discussion has been concerned with monochrome images. Fora colour image the three different colour channels (red, green and blue)may be analysed separately (and corrected separately), using eitherequation (2) or equation (37) set out above.

However an alternative method is preferred, in which a compositefunction is used and a three-variable optimisation process carried out.One such composite function is set out as equation (38):

                                          (38)${S\left( {\lambda_{r},\lambda_{g},\lambda_{b}} \right)} = {\quad{\frac{1}{K}{\left\{ {{\sum\limits_{k = 1}^{k = K}\left( \frac{p_{k_{r}} - \overset{\_}{p_{r}}}{\left( {\overset{\_}{p_{k_{r}}} - \lambda_{r}} \right)} \right)^{2}} + {\sum\limits_{k = 1}^{k = K}\left( \frac{p_{k_{g}} - \overset{\_}{p_{k_{g}}}}{\left( {\overset{\_}{p_{k_{g}}} - \lambda_{g}} \right)} \right)^{2}} + {\sum\limits_{k = 1}^{k = K}\left( \frac{p_{k_{b}} - \overset{\_}{p_{k_{b}}}}{\left( {\overset{\_}{p_{k_{b}}} - \lambda_{b}} \right)} \right)^{2}}} \right\} \cdot \exp}{\left\{ {\frac{1}{K}\frac{\left( {{\sum\limits_{k = 1}^{k = K}{\ln \left( {\overset{\_}{p_{k_{r}}} - \lambda_{r}} \right)}^{2}} + {\sum\limits_{k = 1}^{k = K}{\ln \left( {\overset{\_}{p_{k_{g}}} - \lambda_{g}} \right)}^{2}} + {\sum\limits_{k = 1}^{k = K}{\ln \left( {\overset{\_}{p_{k_{b}}} - \lambda_{b}} \right)}^{2}}} \right)}{(3)}} \right\}.}}}$

where:

-   -   K is the number of pixels to be processed;    -   p_(k) _(r) , p_(k) _(g) , p_(k) _(b) are the values of pixel k        for the red, green and blue channels respectively;

p_(k) _(r) , p_(k) _(g) , p_(k) _(b) are the values of pixel k for thered, green and blue channels respectively after application of thelowpass filter described above;

-   -   (λ_(r), λ_(g), λ_(b)) are the values of the parameters to be        determined; and    -   S(λ_(r), λ_(g), λ_(b)) is the function to be optimised.

As before, subject to certain assumptions, (that the fractionalvariation in brightness is approximately independent of the actualbrightness, and that the fractional variation is similar in the red,green and blue colour channels) the function S(λ_(r), λ_(g), λ_(b)) hasa stationary point (a minimum) when (λ_(r), λ_(g), λ_(b)) indicates thecomponent of each channel attributable to airlight.

The evaluation technique described above can be used for colour imagesusing equation (38) above.

Here, the parameter RS referred to above is replaced by three RSparameters, one for each channel of the image. These are definedaccording to equations (39), (40) and (41) set out below:

$\begin{matrix}{{{{RS}_{r}(i)} = {\frac{1}{M}{\sum\limits_{m = 1}^{M}\frac{\left( {P_{r_{im}} - {\overset{\sim}{P}}_{r_{im}}} \right)^{2}}{\left( {{\overset{\sim}{P}}_{r_{im}} - \lambda_{r_{i}}} \right)^{2}}}}},} & (39) \\{{{{RS}_{g}(i)} = {\frac{1}{M}{\sum\limits_{m = 1}^{M}\frac{\left( {P_{g_{im}} - {\overset{\sim}{P}}_{g_{im}}} \right)^{2}}{\left( {{\overset{\sim}{P}}_{g_{im}} - \lambda_{g_{i}}} \right)^{2}}}}},} & (40) \\{{{{RS}_{b}(i)} = {\frac{1}{M}{\sum\limits_{m = 1}^{M}\frac{\left( {P_{b_{im}} - {\overset{\sim}{P}}_{b_{im}}} \right)^{2}}{\left( {{\overset{\sim}{P}}_{b_{im}} - \lambda_{b_{i}}} \right)^{2}}}}},} & (41)\end{matrix}$

Similarly, the parameter RL referred to above is replaced by threeparameters, again one for each channel of the image, as shown inequations (42), (43) and (44):

$\begin{matrix}{{{{RL}_{r}(i)} = {\frac{1}{3M}{\sum\limits_{m = 1}^{M}{\ln \left( {{\overset{\sim}{P}}_{r_{im}} - \lambda_{r_{i}}} \right)}^{2}}}},} & (42) \\{{{{RL}_{g}(i)} = {\frac{1}{3M}{\overset{M}{\sum\limits_{m = 1}}{\ln \left( {{\overset{\sim}{P}}_{g_{im}} - \lambda_{g_{i}}} \right)}^{2}}}},} & (43) \\{{{{RL}_{b}(i)} = {\frac{1}{3M}{\sum\limits_{m = 1}^{M}{\ln \left( {{\overset{\sim}{P}}_{b_{im}} - \lambda_{b_{i}}} \right)}^{2}}}},} & (44)\end{matrix}$

Similarly again, the parameter RTL is replaced by three values, one foreach channel according to equation (45), (46) and (47) as set out below.

$\begin{matrix}{{{{RTL}_{r}(n)} = {\frac{1}{n}{\sum\limits_{i = 1}^{n}{{RL}_{r}(i)}}}},} & (45) \\{{{{RTL}_{g}(n)} = {\frac{1}{n}{\sum\limits_{i = 1}^{n}{{RL}_{g}(i)}}}},} & (46) \\{{{{RTL}_{b}(n)} = {\frac{1}{n}{\sum\limits_{i = 1}^{n}{{RL}_{b}(i)}}}},} & (47)\end{matrix}$

Processing is carried out for a single channel at a time and accordinglya parameter S is again replaced by three parameters according toequations (48), (49) and (50):

$\begin{matrix}{{{S_{r}(n)} = {{\frac{1}{n}\left\lbrack {\sum\limits_{i = 1}^{n}{{RS}_{r}(i)}} \right\rbrack} \cdot {\exp \left\lbrack {{RTL}_{r}(n)} \right\rbrack}}},} & (48) \\{{{S_{g}(n)} = {{\frac{1}{n}\left\lbrack {\sum\limits_{i = 1}^{n}{{RS}_{g}(i)}} \right\rbrack} \cdot {\exp \left\lbrack {{RTL}_{g}(n)} \right\rbrack}}},} & (49) \\{{S_{b}(n)} = {{\frac{1}{n}\left\lbrack {\sum\limits_{i = 1}^{n}{{RS}_{b}(i)}} \right\rbrack} \cdot {{\exp \left\lbrack {{RTL}_{b}(n)} \right\rbrack}.}}} & (50)\end{matrix}$

Accordingly, equation 37 can then be rewritten as:

S(λ_(R1),λ_(G1),λ_(B1), . . . ,λ_(Rn),λ_(Gn),λ_(Bn))=S_(r)(n)·Exp[RTL_(g)(n)+RTL _(b)(n)]+

S _(g)(n)·exp[RTL_(b)(n)+RTL _(r)(n)]+

S _(b)(n)·exp[RTL_(r)(n)+RTL _(g)(n)]  (51)

And it will be appreciated that values of the variables S and RTL can becomputed as in the single channel case described above. In acomputational implementation the variable RowProdSum is replaced by athree element array and the variable Row_Sum_Log is again replaced by athree element array. The RTL variable RTB is again replaced by anappropriate three element array.

The component of airtight within an image may vary in dependence uponthe distance between the camera, and a point in the scene represented bya pixel. FIG. 4 illustrates a camera 3 positioned to generate an imageof a scene. The camera 3 is directed at an angle θ to the horizontal,and has a field of view angle FOV.

It is known that airtight (λ) will typically vary in dependence upon thedistance between the camera 3 and the image plane 4, referred to as D,in accordance with equation (52):

λ∝(1−exp(−β_(sc)(ω)·D)  (52)

where β_(sc) is the total volume scattering coefficient for light ofwavelength cu. A normalized volume scattering β_(sc)′(ω) coefficient isdefined, such that:

β_(sc)′(ω)=β_(sc)(ω)·q,  (53)

Where q is the value of the distance D corresponding to the centre ofthe image. In the described embodiment of the present invention thenormalized scattering coefficient β_(sc)′ is represented by three scalarvariables X₀, X₁ and X₂, corresponding to nominal wavelengths in thered, green and blue colour bands respectively.

The geometry of the scene illustrated in FIG. 4 is, in general, unknown.Two further scalar variables, X₃ and X₄, are therefore defined asfollows:

$\begin{matrix}{X_{4} = {F\; O\; V}} & (54) \\{X_{3} = \frac{F\; O\; V}{\theta}} & (55)\end{matrix}$

That is:

Where FOV is the field of view angle and θ is the angle between thehorizontal and the distance from the camera to the center of the image.

$\begin{matrix}{\theta = \frac{X_{4}}{X_{3}}} & (56)\end{matrix}$

All five parameters (X₀, X₁ X₂ X₃ X₄) are determined by multi-variableoptimisation.

The distance between the camera and the image plane varies in dependenceupon position within the image plane. This position is represented by avariable ν which takes a value of 0.5 at one extreme, −0.5 at anotherextreme and 0 at the centre.

Applying the sin rule to FIG. 4.

$\begin{matrix}{\frac{D}{\sin \; \theta} = \frac{q}{\sin \; \gamma}} & (57)\end{matrix}$

It can be seen that

γ=180−(θ+φ)  (58)

Given that:

φ=−νX ₄  (59)

then:

γ=180−(θ−νX ₄)  (60)

and:

sin γ=sin(θ−νX ₄)  (61)

So

$\begin{matrix}{D = {\sin \; \theta \frac{q}{\sin \; \gamma}}} & (62) \\{D = {q\frac{\sin \; \theta}{\sin \left( {\theta - {v\; X_{4}}} \right)}}} & (63) \\{D = {q{\frac{\sin \frac{X_{4}}{X_{3}}}{\sin \; {X_{4}\left( {\frac{1}{X_{3}} - V} \right)}}.}}} & (64)\end{matrix}$

Substituting equation (64) into equation (52) gives the airlight of apixel for a monochrome image or a single channel (e.g. the red channel)of a colour image;

$\begin{matrix}{{{\lambda \left( {x_{0}x_{3}x_{4}} \right)} = {C\left\{ {1 - {\exp\left\lbrack {{- x_{0}}\frac{\sin \left( {x_{4}/x_{3}} \right)}{\sin \left( {x_{4}\left\lbrack {\frac{1}{x_{3}} - v} \right\rbrack} \right)}} \right\rbrack}} \right\}}},} & (65)\end{matrix}$

where C is the constant of proportionality in equation (52). The termλ(X₀, X₃, X₄) then replaces λ in equation (2) or (37). The value ofconstant C depends on the product of the atmospheric radiance and thecamera gain. In general these parameters are not known. An estimate forC may be formed in various ways. One way is to use the highest pixelvalue in the image.

If the colour version of the equation (equation (38)) is used then X₀ isreplaced with X₁ or X₂ for the green and blue channels as appropriate.This term is computed separately for each row of pixels in the image.

When the field of view angle is relatively small, approximately 15degrees or less, the following approximation may be made.

$\begin{matrix}{{{\lambda \left( {x_{0}x_{3}x_{4}} \right)} \approx {\lambda \left( {x_{0}x_{3}} \right)}} = {C{\left\{ {1 - {\exp\left\lbrack {{- x_{0}}\frac{1/x_{3}}{\left\lbrack {\frac{1}{x_{3}} - v} \right\rbrack}} \right\rbrack}} \right\}.}}} & (66)\end{matrix}$

This expression involves only one geometrical parameter. This shows thatthe field of view angle (parameter X₄ in equation (67)) has littleeffect in practice. In the preferred implementation equation (56) isused but with a fixed value for X₄. If the field of view of the camerais not known, then X₄ is set to some small value, say 1 degree.

The preceding discussion has been concerned with estimating airlight instatic images. It will be appreciated that the techniques describedabove can be applied to video data comprising a plurality of frames, byprocessing each frame in turn. When video data is processed in this waycare must be taken to ensure that the image processing operations aresufficiently fast so as not to cause artefacts of “jerkiness” in theprocessed video data. Thus, it may be necessary to compromise imageprocessing quality for the sake of speed.

In accordance with some embodiments of the present invention imageprocessing is carried out using an arrangement illustrated in FIG. 5.Colour input video data 6 in which each pixel of each frame isrepresented by separate red, green and blue channels is received andpassed concurrently to an analysis module 7 and a video processingmodule 8. The analysis module operates relatively slowly and isconfigured to generate values for the parameters c and m of equation(1). Determination of c may be carried out, for example, using themethod outlined above. Values of the parameters c and in are stored in acoefficients buffer 9.

Not all of the image data may be required for the image analysis. Forexample the image may be sub-sampled by a factor of two such thatanalysis is based upon only half of the pixel value. This saves time inthe computation. Also the analysis is usually carried out for apre-defined part of the image in order to avoid borders or screenannotation (which is sometimes used in CCTV installations).

The video processing module 8 operates by processing each pixel of areceived frame in accordance with equation (1), reading appropriatevalues for the parameters c and m from the coefficients buffer 9. Dataprocessed by the video processing module 8 is output to a display screen10.

Given that the analysis module 7 operates slowly compared to the videoprocessing module 8, the video processing module 8 will typically applycoefficients generated using an earlier frame of video data. Given thatchanges in airlight tend to happen over a prolonged period of time, thisis found in practice to have little effect on image quality.

The analysis module 7 preferably produces values of c and m for eachpixel of the processed frame. That is a method taking into accountvariation of airlight with distance from a camera position such as thatdescribed above or that described in EP0839361 is preferably used. Insuch embodiments of the invention a total of six tables are used, twofor each colour channel. Each table has one entry for each pixel in theimage. The red, green and blue channels are processed separately andthen re-combined to form the final image.

The value of m can be adjusted to provide compensation for unevenillumination in addition to the airlight compensation if required.Several methods to estimate local illumination levels are available inthe literature, such as that described in, for example Jobson, J.,Rahman, Z., Woodell, G. A., “A multiscale retinex for bridging the gapbetween colour images and the human observation of scenes”, IEEETransactions on Image Processing, Vol. 6, Issue 7, July 1997, pp.965-976.

Illumination estimation and compensation is now briefly described.Illumination compensation is based on an idea that illumination changesslowly, and that estimates of illumination distribution can be generatedusing a low pass filter. It is known that low pass filters such as thehomomorphic filter can be used to estimate illumination for the purposesof image enhancement. This is described in, for example, Jobson, Rahmanand Woodell referenced above.

Illumination estimation methods of the type described above are based onthe Lambertian model for diffuse reflection. This is described in AnyaHurlbert, “Formal connections between lightness algorithms”, J. Opt.Soc. Am. A, Vol. 3, No. 10, pp. 1684-1693, October 1986. If the visiblesurfaces in an image have a Lambertian characteristic and theillumination is confined to a narrow range of wavelengths then:

$\begin{matrix}{{I\left( {x,y} \right)} = {{R\left( {x,y} \right)}{L\left( {x,y} \right)}}} & (67) \\{{R\left( {x,y} \right)} = \frac{I\left( {x,y} \right)}{L\left( {x,y} \right)}} & (68)\end{matrix}$

where I(x,y) is the intensity value at pixel (x,y), R(x,y) is the scalarreflectance factor and L(x,y) is a scalar irradiance value. The imageprocessing procedure is simple; first form an estimate for L(x,y) andthen scale by c/L(x,y), where the constant c is chosen to keep thepixels within the display range.

This model can be generalized to colour images under the assumption thatnarrow-band filters are used for each colour channel. Even when thespectrum is localized using conventional ROB filters the model definedby (68) gives an approximation to the true response. The benefits ofsuch decomposition include the possibility of removing illuminationeffects of back/front lighting, enhancing shots that include spatiallyvarying illumination such as images that contain shadows.

Various methods have been proposed for estimation of the illuminationfunction L(x,y). In Akira Suzuki, Akio Shio, Hiroyuki Arai, and SakuichiOhtsuka, “Dynamic shadow compensation of aerial images based on colorand spatial analysis”, Pattern Recognition, 2000. Proceedings. 15thInternational Conference on, 3-7 Sep. 2000, Vol. 1, pp. 317-320, amorphological smoothing filter is used. The homomorphic filter,described by Thomas Stockham, Jr. in “Image processing in the context ofa visual model,”, Proc. IEEE, Vol. 60, No. 7, pp. 828-842, 1972, is anearly illumination compensation technique. The retinex algorithm,described by Jobson, Rahman, and Wooden in “Properties and performanceof a centre/surround retinex”, IEEE Transactions on Image Processing,Vol. 6, No. 3, 1997, pp. 451-462, is a more recent method.

The illumination estimation and compensation techniques described abovecan be applied in the context of a system such as that illustrated inFIG. 5. Specifically, the analysis module 7 can estimate illumination,while the video processing module 8 can process data to performillumination compensation.

Although preferred embodiments of the invention have been describedabove, it will be appreciated that various modifications can be made tothose embodiments without parting from the spirit and scope of thepresent invention as defined by the amended claims.

The description set out above has been concerned with image processingtechniques concerned with removing the effects of airlight from animage. It will be appreciated that there are other sources of imageoffsets apart from the airlight effect. In particular image offsets canbe generated by the thermal response of a camera—the so-called “darkcurrent”. Offsets can also be generated by dirt on the camera lens orprotective window. Additionally noise can be caused by so called“pedestal” errors. Some video standards, including certain variants ofthe PAL standard, use a non-zero voltage level to represent black. Thislevel is known as the “Pedestal”. If the video decoder makes aninappropriate allowance for the pedestal the effect is to shift all ofthe brightness levels either up or down. All of these types of offsetscan be detected (and mitigated) by the methods described above.

1. A method of processing a plurality of frames of video data togenerate enhanced video data, comprising: processing a first frame ofvideo data to generate an estimate of noise included within said frameof video data; wherein said noise is at least partially attributable toatmospheric backscattered light; storing data indicative of saidestimate of noise; processing at least one further frame of video data,which further frame indicates the same scene as said first frame ofvideo data; to remove noise from said further frame of video data bysubtracting said estimate of noise included in said first frame of videodata and; outputting an enhanced frame of video data generated usingsaid further frame of video data.
 2. A method according to claim 1,wherein said subtracting comprises subtracting said estimate from valuesof said pixels of said video frame.
 3. A method according to claim 2,further comprising rescaling values of said pixels after subtractingsaid estimate from values of said pixels.
 4. A method according to claim1, wherein processing at least one further frame of video data comprisesprocessing a plurality of further frames of video data using said storedestimate of noise; and the method comprises outputting a plurality ofenhanced frames of video data.
 5. A method according to claim 1, whereinsaid processing said first frame of video data to generate said estimatecomprises: evaluating a function for different values of said estimate,said function taking as input only an estimate of said noise; anddetermining an estimate of said noise for which said function has anoptimum value.
 6. A method according to claim 5, wherein first frame ofvideo data comprises pixel values for each of a plurality of pixels ofsaid frame, and said function takes as input as least a subset of saidpixel values.
 7. A method according to claim 5, further comprising:filtering said video data to generate filtered video data comprisingfiltered pixel values; wherein said function takes as input at least asubset of said filtered pixel values.
 8. A method according to claim 5,wherein said optimum value of said function is a stationary point ofsaid function.
 9. A method according to claim 8, wherein said stationarypoint is a minimum.
 10. A method according to claim 5, wherein saidfunction is of the form:${S(\lambda)} = {\frac{1}{K}{\sum\limits_{k = 0}^{k < K}{{\left( \frac{p_{k} - \overset{\_}{p_{k}}}{\overset{\_}{p_{k}} - \lambda} \right)^{2} \cdot \exp}\frac{1}{K}{\sum\limits_{k = 0}^{k < K}\left( {\ln \left( {\overset{\_}{p_{k}} - \lambda} \right)}^{2} \right)}}}}$where: K is the number of pixels to be processed; p_(k) is the value ofpixel k p_(k) is the value of pixel k after application of the lowpassfilter described above; λ is the value of the estimate; and S(λ) is thefunction to be optimized.
 11. A method according to claim 5, whereinsaid optimum value is determined using a numerical analysis technique.12. A method according to claim 10, wherein said estimate is initiallyzero.
 13. A method according to claim 5, comprising generating aplurality of estimates, each estimate being an estimate of noise in arespective subset of pixels of said video frame.
 14. A method accordingto claim 1, wherein said step of processing at least one further framecomprises processing each pixel of said further video frame to removesaid estimate of said noise to generate output pixel values.
 15. Amethod according to claim 14, further comprising: multiplying saidoutput pixel values by a predetermined coefficient.
 16. A methodaccording to claim 1, wherein said noise is at least partiallyattributable to dirt on a camera lens.
 17. A method according to claim1, wherein said noise comprises an offset applied to at least a subsetof pixels of said video frame, said offset being substantially equal forall pixels of said subset pixels of said video frame.
 18. A methodaccording to claim 17, wherein said offset is applied to all pixels ofsaid video frame.
 19. A method according to claim 1, wherein said noisecomprises an offset applied to at least some pixels of said video frame,said offset being determined for particular independence upon a distancebetween a camera position used to generate said video frame, a pointrepresented by that pixel.
 20. A method according to claim 1, whereinsaid processing said first frame of video data further generates acontrast enhancement parameter.
 21. A computer apparatus for processinga plurality of frames of video data comprising: a program memorycontaining processor readable instructions and; a processor configuredto read and execute instructions stored in said program memory: whereinsaid processor readable instructions comprise instructions configured tocause the computer to carry out a method according to claim 1.